#!/usr/bin/perl
use strict;
use warnings;
use Getopt::Long;
use Cwd;

my $list;
my @all_sample;
my @sample_name;

my $force=0;
GetOptions(
'list=s' => \$list,

'force!' => \$force,
);


my $usage=<<USAGE;
-list          <str>       Sample_list
USAGE

if (!$list) {
	print $usage;
	exit;
}

my $pwd_path = getcwd();

delete_windows_enter($list);


####main####
makefolder();
###1. 检查文件是否齐全
	my @report_back = @{check()};
	my %report_hash = %{$report_back[0]}; ###拿到HASH
	my $report_sample = $report_back[1]; ###拿到实例用的sample

###2. 整合文档并转化格式
	get_table(); ###将table合并起来并转化格式


###3. 出报告
	get_png();
	report();




####sub#####
sub report{
	my $report=<<REPORT;
![Alt text](./logo.png)

#美吉生物结题报告 

> **人外显子测序及分析项目**

----------
<br/>
 **客户：xx**            
 **2016 年xx 月xx 日**
<br/>
 ----------
<pre><code>
美吉生物结题报告
一. 工作流程
	1.1 测序实验流程
	1.2 生物信息分析流程
二. 项目结果报告
	2.1 原始测序数据统计
	2.2 原始测序数据质控
	2.3 比对结果数据统计
	2.4 个体SNP和InDel检测
	2.5 突变率和突变类型分析
	2.6 个体SNP和InDel的注释
三. 附录
	3.1 结果文件列表
	3.2 文件解压缩方法
	3.3 文件打开或浏览方法
	3.4 参考文献
	3.5 联系方式
<\/code><\/pre>
----------
##一. 工作流程
###1.1 测序实验流程
每个样本的基因组DNA随机打断成150~200bp的片段并构建文库，参照SureSelectXT Target Enrichment System for Illumina Paired-End Sequencing Library, Illumina HiSeq and MiSeq Multiplexed Sequencing Platforms, Protocol Version 1.3.1 进行样本捕获及制备，过程如图1-1所示：
<center>![Alt text](./Agilent.png)

图1-1  Agilent 外显子捕获流程</center>

###1.2 生物信息分析流程
生物信息分析流程步骤如下：
<br>**1.** 用clean data与参考基因组Hg19比对;
<br>**2.** 利用Picard-tools去除PCR-duplication产生的测序reads;
<br>**3.** 根据比对结果，计算exome区每个位点的覆盖深度、总体平均深度和覆盖的百分比，同时计算捕获效率；
<br>**4.** 利用GATK软件对InDel附近的reads进行局部的重新对齐，获得realign后的bam文件，可以消除InDel周边的假阳性SNP；
<br>**5.** 使用校正后的bam文件，结合Agilent Exon V5+UTR信息，利用GATK软件对exome区进行SNP、small InDel的检测，再根据项目的研究目的，结合人类dbsnp信息，过滤掉人群SNP位点（通常与疾病无关）。
<br>**6.** 利用Annovar软件[2]和参考基因组hg19 refGene对突变进行注释，得到SNP、small InDel的注释信息。接下来利用SIFT、Ployphon-2等软件对筛选出的潜在突变进行氨基酸置换的预测，以判断突变对蛋白质结构和功能的影响。在分析一些复杂疾病时，才利用PLINK软件计算突变在疾病样本和正常样本中次等位基因的频率，用Fisher Exact Test 判断突变位点与疾病是否存在相关性。

生物信息分析流程如图1-2：



<center>![Alt text](./analysis_pipeline_1.png)


图1-2  生物信息分析流程</center>

##二.  项目结果报告
###2.1 原始测序数据统计

Illumina Hiseq 测序得到的原始图像数据经过Base Calling转化为序列数据，以FASTQ文件格式来存储，是用户得到的最原始的数据文件。FASTQ格式文件包含reads的序列信息以及reads的测序质量信息。测序数据随机截取结果如下所示：
<br>\@HWI-ST531R:144:D11RDACXX:4:1101:1212:1946 1:N:0:ATTCCT
<br>ATNATGACTCAAGCGCTTCCTCAGTTTAATGAAGCTAACTTCAATGCTGAGATCGTTGA
<br>+HWI-ST531R:144:D11RDACXX:4:1101:1212:1946 1:N:0:ATTCCT
<br>?A#AFFDFFHGFFHJJGIJJJIICHIIIIJJGGHIIJJIIJIIJIHGI\@FEHIIJBFFH

每条read包含4 行信息，其中第一行和第三行由文件识别标志和读段名（ID）组成（第一行以“@” 开头而第三行以“+” 开头；第三行中 ID可以省略，但“+” 不能省略），第二行为碱基序列，而第四行是第二行中的序列内容每个碱基所对应的测序质量值。

得到原始测序数据后，对测序结果进行质量统计，原始数据的质量统计结果见表2-1。

表2-1  原始测序数据质量统计表

$report_hash{"raw_Q20"}

###2.2 原始测序数据质控

Illumina 测序属于第二代测序技术，单次运行能产生数十亿级的reads，如此的海量数据无法逐个展示每条read的质量情况；生物信息运用统计学的方法，对所有测序reads进行碱基分布和质量波动统计，可以从宏观上直观地反映出测序样本的质量和文库构建质量。
根据测序数据统计了外显子文库的片段大小，一般外显子文库片段大小在100到300碱基之间。
<br>图2-1为测序片段长度分布图；
<br>图2-2为原始数据碱基组成分布图(又称为GC偏差图)；
<br>图2-3为原始数据碱基质量分布图；
<br>图2-4  reads的quality的均值分布；
<br>图2-5  reads的duplication的比例。
<br>结果目录为**fastqc**，目录中文件说明如下：其中包括测序数据统计表，各样本的fastqc（样本名/*.fastqc）的统计结果。

<center><img src="./$report_hash{'insert_plot'}" width = "700" height = "400">

图2-1 测序片段长度分布图</center>

<br/>
<center><img src="./$report_hash{'ATCG_per'}" width = "700" height = "400">

图2-2  原始数据碱基组成分布例图

**注：**横坐标是reads碱基坐标，纵坐标是所有reads的A、C、G、T、N碱基分别占的百分比。每个位置上，A、C、G、T在开始有所波动，后面会趋于稳定。一般情况下A与T相等，C与G相等，各碱基所占百分比会因物种差异而不同。基因组项目中，建库比较均匀的情况下，代表不同碱基的四种颜色的分界线应该波动较小，比例较为一致。</center>

从上图可知，该文库碱基分布均匀，可用于后续分析。
<br/>
<center><img src="./$report_hash{'base_quality_per'}" width = "700" height = "400">

图2-3  原始数据碱基质量分布例图</center>

**注：**横坐标是reads碱基坐标，纵坐标是reads的碱基质量，图中垂直线指定的范围是所有reads碱基的综合质量，黄色垂直方块是质量的四分位值范围。

从上图可知，该文库碱基质量良好，可用于后续分析。

<center><img src="./$report_hash{'seq_quality_per'}" width = "700" height = "400">

图2-4  reads的quality的均值的分布

**注：**每条reads的quality的均值的分布，横轴为quality，纵轴是reads数目。</center>

从上图可知，该文库reads碱基质量良好，可用于后续分析。

<center><img src="./$report_hash{'reads_dup'}" width = "700" height = "400">

图2-5  reads的duplication的比例

**注：** duplication指的是测序完全一样的reads，图2-5统计序列完全一样的reads的频率。横坐标是duplication的次数，纵坐标是duplicated reads的数目，以unique reads的总数作为100%。测序深度越高，越容易产生一定程度的duplication，这是正常的现象，但如果duplication的程度很高，就提示我们可能有bias的存在（如建库过程中的PCR duplication）。上图的情况中，duplication的比例<10%，该文库数据质量质量良好，可用于后续分析。 </center>

###2.3 比对结果数据统计


采用Illumina Hiseq测序技术对样品的DNA进行paired-end(PE)测序，由于Illumina Hiseq的原始测序数据会存在一些质量比较低的数据，为了使后续的组装更加准确，会对原始数据进行质量剪切，具体步骤如下：
<br>1) 去除reads中的adapter序列； 
<br>2) 剪切前去除两端含有非AGCT的碱基；
<br>3) 修剪测序质量较低的reads末端 (测序质量值小于Q20)；
<br>4) 去除N的数目大于等于5个（默认）的reads；
<br>5) 舍弃修剪后长度小于25bp的小片段。

质量剪切后得到clean data，并对clean data进行质量统计，clean data的质量统计结果见表2-2。

表2-2  clean data测序数据质量统计表

$report_hash{"clean_Q20"}


使用BWA软件将clean data的片段比对回参考基因组，并使用Picard软件去除PCR-duplication产生的测序序列。本项目中采用的参考基因组是hg19。人的基因组大小约3G，Agilent捕获外显子区域大小约为75Mb。
比对结果的统计见表2-3和表2-4：

表2-3  参考基因组比对结果数据统计表

$report_hash{"mapping_statistic"}

<br/>
表2-4  外显子区域比对结果数据统计表

$report_hash{"exon_statistic"}

###2.4 个体SNP和InDel检测

在SNP和InDel calling之前，先对比对的bam文件进行处理：
InDel附近的alignment通常有一定误差，需要利用已知的InDel信息进行realignment来对 InDel附近的reads进行局部的重新对齐，获得realign后的bam文件，可以消除indel周边的假阳性SNP。最后再结合dbsnp信息和Agilent Exon V5+UTR的外显子区间信息进行SNP和InDel的检测。

使用的dbsnp信息如下：
<br> ##dbSNP_BUILD_ID=138
<br> ##fileDate=20130806
<br> ##phasing=partial
<br> ##reference=file:///humgen/gsa-hpprojects/GATK/data/ucsc.hg19/ucsc.hg19.fasta
<br> ##source=dbSNP
<br> ##variationPropertyDocumentationUrl=ftp://ftp.ncbi.nlm.nih.gov/snp/specs/dbSNP_BitField_latest.pdf 

在个体SNP和InDel的检测过程中，要求该位点的测序深度大于等于5（默认）。

结果文件为**annot_result_final.xls_cancer_gene.xls**
vcf文件格式可直接用记事本或Notepad++等编辑器打开，具体解释如下表2-5。

表2-5 vcf文件各列的含义

| 简称 | 具体解释|
| :--------: |:--------:|
|CHOM POS|参考序列名和variant位置，若为InDel，位置为InDel的第一个碱基位置|
|ID|Variant的ID，比如在dbSNP中有该SNP的id，则会在此行给出；若没有，则用’.'表示其为一个novel variant|
|REF ALE|参考序列的碱基和Variant的碱基|
|QUAL|Phred格式(Phred_scaled)的质量值，表示在该位点存在variant的可能性；该值越高，则variant的可能性越大；计算方法：Phred值 = -10 * log (1-p) p为variant存在的概率; 通过计算公式可以看出值为10的表示错误概率为0.1，该位点为variant的概率为90%|
|Gene|注释的官方基因名称|
|transcript|转录本名称|
|Fun|对变异位点所在的区域进行注释（exonic, splicing, UTR5, UTR3, intronic, ncRNA_exonic, ncRNA_intronic, ncRNA_UTR3, ncRNA_UTR5, ncRNA_splicing, upstream, downstream, intergenic）|
|AAChange|氨基酸改变，只有当Func列为exonic或exonic;splicing时，该列才有结果。按照每个转录本进行注释（例如，NADK:NM_001198995:exon10:c.1240_1241insAGG:p.G414delinsEG，其中，NADK表示该变异所在的基因名称，NM_001198995表示该变异所在的转录本ID，exon10表示该变异位于转录本的第10个外显子上，c.1240_1241insAGG表示该变异引起cDNA在第1240和1241位之间插入AGG，p.G414delinsEG表示该变异引起蛋白序列在第414位上的氨基酸由Gly变为Glu-Gly；FMN2:NM_020066:exon1:c.160_162del:p.54_54del，表示该变异引起cDNA的第160到162位发生删除，p.54_54del表示该变异引起蛋白序列在第54位上的氨基酸删除）|
|GeneDetail|描述UTR、splicing、ncRNA_splicing或intergenic区域的变异情况。当Func列的值为exonic、ncRNA_exonic、intronic、ncRNA_intronic、upstream、downstream、upstream;downstream、ncRNA_UTR3、ncRNA_UTR5时，该列为空；当Func列的值为intergenic时，该列格式为dist=1366;dist=22344，表示该变异位点距离两侧基因的距离|
|ExonicFunc|外显子区的SNV or InDel变异类型（SNV的变异类型包括synonymous_SNV, missense_SNV, stopgain_SNV, stopgloss_SNV和unknown；Indel的变异类型包括frameshift insertion, frameshift deletion, stopgain, stoploss, nonframeshift insertion, nonframeshift deletion和unknown）|
|cytoBand|该变异位点所处的染色体区段（利用Giemas染色观察得到的）。如果变异位点跨过多个区段，用短横线连接|
|phyloP46way|给出该变异位点是否位于保守区域中；如果在保守区域中，给出该区域的分值。保守区域是由phastCons程序基于脊椎动物全基因组比对预测得到的，46way是指使用的物种数目为46个。有两个值，一是score，是该保守区域的分值（transformed log-odds scores），取值为0-1000；二是name，是ANNOVAR给该保守区域的名称|
|SuperDups|检测该变异位点是否位于重复片段（segmental duplication）中。重复区域中检测到的遗传变异大多数是由于序列比对错误造成的，所以被注释到segmental duplications的变异需要谨慎对待，很可能是假阳性位点。给出两个值，一是Name，表示基因组中与该变异位点所在区域相似的片段的位置；二是Score，表示两个相似片段的序列一致性。例如，Score=0.994828;Name=chr19:60000，表示chr19:60000所在片段跟该变异位点所在片段相似，序列一致性为0.994828|
|esp6500siv2|国家心肺和血液研究所外显子组测序计划（NHLBI-ESP project，esp6500si_all数据库中包含SNP变异、Indel变异和Y染色体上的变异）的所有个体中，突变碱基的等位基因频率（alternative allele frequency）|
|1000g2014oct_all|给出千人基因组计划数据（2014年10月公布的版本）的所有人群中，该变异位点上突变碱基的等位基因频率|
|1000g2014oct_afr|给出千人基因组计划数据（2014年10月公布的版本）的非洲人群中，该变异位点上突变碱基的等位基因频率|
|1000g2014oct_eas|给出千人基因组计划数据（2014年10月公布的版本）的亚洲人群中，该变异位点上突变碱基的等位基因频率|
|1000g2014oct_eur|给出千人基因组计划数据（2014年10月公布的版本）的欧洲人群中，该变异位点上突变碱基的等位基因频率|
|snp138|该变异在dbSNP（版本138）中的ID|
|Genotype|基因型|
|INFO/INDEL|INDEL标示，如果是1则为INDEL|
|DEPTH|深度信息，分别是REF/ALT的深度信息，用逗号隔开|
|ljb26_sift|SIFT分值，表示该变异对蛋白序列的影响，包含两个值，一是SIFT初始分值，二是T或者D。当该变异同时影响多个蛋白序列时，对每条蛋白序列有一个SIFT值，取最小值。SIFT分值越小越“有害”，表明该SNP导致蛋白结构或功能改变的可能性大；D: Deleterious (sift<=0.05); T: tolerated (sift>0.05)）|
|ljb26_hdiv|利用PolyPhen2基于HumanDiv数据库预测该变异对蛋白序列的影响，用于复杂疾病。该列包含两个值，第一个是PolyPhen 2分值，数值越大越“有害”，表明该SNP导致蛋白结构或功能改变的可能性大；第二个是D或P或B（D: Probably damaging (>=0.957), P: possibly damaging (0.453<=pp2_hdiv<=0.956); B: benign (pp2_hdiv<=0.452)）|
|ljb26_hvar|利用PolyPhen2基于HumanVar数据库预测该变异对蛋白序列的影响，用于单基因遗传病。该列包含两个值，第一个是PolyPhen 2分值，数值越大越“有害”，表明该SNP导致蛋白结构或功能改变的可能性大；第二个是D或P或B（D: Probably damaging (>=0.909), P: possibly damaging (0.447<=pp2_hvar<=0.909); B: benign (pp2_hvar<=0.446)）|
|ljb26_lrt|LRT分值，表示该变异对蛋白序列的影响，包含三个值，一是LRT初始分值，二是D、N或者U(D: Deleterious; N: Neutral; U: Unknown)。第一个值越大越“有害”，表明该SNP导致蛋白结构或功能改变的可能性大|
|ljb26_mt|MetaLR分值，表示该变异对蛋白序列的影响，包含两个值，一是MetaLR初始分值，二是D或者T(D: Deleterious; T: Tolerated)|
|cancer_gene|基于TCGA、COSMIC数据库的主要肿瘤突变基因|
|HGMD_OMID|人类基因突变数据库_在线人类孟德尔遗传数据库|
|Child_inherited_disease[4]|448种侧重于严重、隐性的儿科疾病|

<br/>
该部分的统计结果见表2-6； SNP与InDel在各染色体目标区域的概率图见图2-6；InDel中Insertion 和Deletion的长度分布统计图见图2-7。

表2-6  个体SNP和InDel结果统计表

$report_hash{"snp_indel_statistic"}
<br/>
<center>
<figure class="half">
 <img src="./$report_hash{'snp_base'}" width = "440" height = "440">
 <img src="./$report_hash{'indel_base'}" width = "440" height = "440">
</figure>

图2-6  SNP与InDel在各染色体目标区域的概率图</center>

<br/>
<center><img src="./$report_hash{'int_del_plot'}" width = "700" height = "400">

图2-7  InDel中Insertion 和Deletion的长度分布统计图</center>

###2.5 突变率和突变类型分析

碱基替换类型可以作为发生在疾病中的特定突变机制；一些疾病中，某类碱基突变类型占绝大多数。所以我们对各个SNP位点发生的转换颠换分布进行了统计，如表4-11：

表2-7 转换颠换类型分布

$report_hash{"mutation_statistic"}

###2.6 个体SNP和InDel的注释

采用Annovar程序结合UCSC Genome Browser的refGene注释信息，对SNP和InDel变异进行注释。对于在exonic区域的SNP和InDel突变位点，将会对突变位点对蛋白翻译的影响进行注释。

结果文件为**annot_result_final.xls_cancer_gene.xls**:

该部分的统计结果见表2-8至2-11：

表2-8  个体SNP注释结果统计表

$report_hash{"snp_annot_statistic"}

<center><img src="./$report_hash{'snp_annot_pie'}" width = "600" height = "600">

图2-8  SNP在不同基因区域的分布统计图

<img src="./$report_hash{'indel_annot_pie'}" width = "600" height = "600">

图2-9  INDEL在不同基因区域的分布统计图</center>

表2-9  个体InDel注释结果统计表

$report_hash{"indel_annot_statistic"}
**表2-8和 表2-9 具体描述和说明参考（http://www.openbioinformatics.org/annovar/annovar_gene.html）

表2-10  exome区域的个体SNP对蛋白翻译影响结果统计表

$report_hash{"snp_pro_statistic"}


表2-11  exonic区域的个体InDel对蛋白翻译影响结果统计表

$report_hash{"indel_pro_statistic"}
  **表2-10 和 表2-11 中unknown是由于注释文件中基因的ORF不完整导致无法进行判断


###2.7 个体SNP对蛋白结构及功能影响注释

通常，在编码区内的SNP位点常与个体之间的表型差异有关。特别是非同义的SNP位点（Nonsynonymous SNP）会导致编码的氨基酸产生改变。对多序列比对结果和3D结构数据的分析，可以得到氨基酸突变对蛋白质结构和功能的影响。

SIFT （Sorting Intolerant From Tolerant，http://sift.jcvi.org/www/SIFT_help.html#SIFT）是一个可以预测氨基酸突变对蛋白质影响的软件。通过目标序列与其近似序列的比对（PSI-BLAST）可以知道目标序列中氨基酸残基的保守程度，SIFT就是通过氨基酸的保守程度进行预测的。

PolyPhen2(http://genetics.bwh.harvard.edu/pph2/dokuwiki/overview)也是一个预测氨基酸突变对蛋白质结构和功能影响的软件。该软件进行预测参考的信息：1、序列数据，依靠序列数据确定突变发生的特异位点，并根据位点的特异性判断对分子结构和功能的影响；2、PSIC值，用来衡量突变在蛋白质家族的罕见程度；3、分子二级结构及相关信息，包括二级结构、可溶性表面积及Phi-psi二面角；4、蛋白质3D结构，将突变氨基酸替换到蛋白质分子的3D结构中看该突变是否会破坏蛋白质的疏水核心。

表2-12 SIFT预测结果的统计表

$report_hash{"sift_statistic"}

**SIFT的预测： 
“DAMAGING”： 突变对蛋白质结构和功能具有一定的破坏性； 
“TOLERATED”： 突变对蛋白质结构和功能影响较小，在可承受范围内；

表2-13 PolyPhen2预测结果的统计表

$report_hash{"hdiv_statistic"}

**PolyPhen2的预测分为三个等级：“benign”，“possibly damaging”，“probably damaging”，危害程度递增。


##三. 附录
###3.1 文件解压缩方法
所有提供的文件均为Linux系统下的文件，压缩包使用“tar -zcvf ”命令压缩，以下为不同系统用户解压缩的方法：
Unix/Linux/Mac用户：使用tar -zcvf *.tar.gz命令
Windows用户：使用WinRAR软件解压缩
###3.2 文件打开或浏览方法
如果在本附录中无特殊说明，所有提供的文件均为Linux系统下文本文件，Unix/Linux用户可以使用more或less命令查看文本
文件内容。对于Windows用户，一般文本文件可以使用写字板或excel打开。推荐使用开源文本编辑器gedit for win32版本
(http://projects.gnome.org/gedit/)或者商业文本编辑器UltraEdit。当文件比较大时，打开文件可能导致Windows系统死机，
建议使用性能较好的计算机或者使用更适合处理大量数据的Unix/Linux系统打开。
数据中可能包含部分图像文件，一般图像文件后缀名为.png、.jpg、.gif等，对于图像文件，Windows用户可以使用图片浏览
器打开，Linux/Unix用户使用display命令打开。
后缀名为svg的文件为文本格式描述的图像文件，Windows用户需要安装Adobe Illustrator软件打开。Linux/Unix用户可以使
用rsvg-view命令查看。公司默认提供“pdf“格式的矢量图，可利用”Adobe Illustrator”软件对该格式图片进行编辑。
Linux下的表格均为制表符(Tab)分割的文本，为了便于阅读，建议使用excel或openoffice等办公软件用表格形式打开，打开
时请用“Tab分割”方式。
###3.3 参考文献
<br>[1] Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics, 1754-1760.2009
<br>[2] Wang K, Li M, Hakonarson H. ANNOVAR: Functional annotation of genetic variants from next-generation sequencing data Nucleic Acids Research, 38:e164, 2010
<br>[3] Carol Jean Saunders et al.: Rapid Whole-Genome Sequencing for Genetic Disease Diagnosis in Neonatal Intensive Care Units Sci Transl Med 4, 154ra135 (2012);
###3.4 联系方式
上海总部
地址： 上海市浦东新区康新公路3399号时代医创园3号楼
电话： 021-51875086
E-mail： meddna\@majorbio.com<br/>
广州分公司
地址： 广州市海珠区荔福路68号广州市微生物所六楼西
电话： 020-61130189
E-mail： seqgz\@majorbio.com<br/>
北京分公司
地址： 北京市海淀区安宁庄东路18号光华创新园新科研楼6楼
电话： 010-51293026,010-51293126
E-mail： seqbj\@majorbio.com
REPORT

	open OUT_SUB,">5_report/report.md" || die $!;
	print OUT_SUB $report;
	close OUT_SUB;
}




sub get_png{
	###子程序的作用是将PNG提取到report下，然后可以很方便的插入结果并体现在markdown下
	foreach my $keys (keys %report_hash) {
		if ($report_hash{$keys} =~ /png$/) {
			system("cp $report_sample/$report_hash{$keys} 5_report");
			$report_hash{$keys} = (split(/\//,$report_hash{$keys}))[-1];
		}
	}
	system("cp /mnt/ilustre/users/kefei.huang/program/exon_pipeline/report_file/logo.png /mnt/ilustre/users/kefei.huang/program/exon_pipeline/report_file/Agilent.png /mnt/ilustre/users/kefei.huang/program/exon_pipeline/report_file/analysis_pipeline_1.png 5_report/");
}



sub get_table{
	foreach my $keys (keys %report_hash) {
		if ($report_hash{$keys} =~ /xls$/) {		
			if ($keys eq "snp_annot_statistic" || $keys eq "indel_annot_statistic" || $keys eq "snp_pro_statistic" || $keys eq "indel_pro_statistic" || $keys eq "sift_statistic" || $keys eq "hdiv_statistic") {
				func_merge($keys);
			}
			if ($keys eq "mapping_statistic" || $keys eq "exon_statistic" || $keys eq "snp_indel_statistic" || $keys eq "mutation_statistic" || $keys eq "clean_Q20" || $keys eq "raw_Q20") {
				row_merge($keys);
			}
		}
	}
}
	sub func_merge{
		my $type = shift;
		my %hash_func;
		@{$hash_func{"header"}} = "type";
		my $count = 0;
		foreach my $samples (@all_sample) {
			push @{$hash_func{"header"}},$sample_name[$count];
			open EACH,"$samples/$report_hash{$type}" || die $!;
			while (<EACH>) {
				chomp;
				my @lines = split(/\t/);
				${$hash_func{$lines[0]}}[$count]=$lines[1];
			}
			$count ++;
		}

		my @header = @{$hash_func{"header"}};
		delete $hash_func{"header"};
		###1.添加header
			my $header = join("|",@header);
			$header = "|$header|";
			$header =~ s/\n\|/\|\n/g;
		###2.添加format
			my $format = "|";
			foreach my $run (@header) {
				$format = "$format:-:|";
			}
		###3.添加其他的东西
			my $content = "$header\n$format\n";
			foreach my $keys (sort keys %hash_func) {
				my @out;
				push @out,$keys;
				push @out,@{$hash_func{$keys}};
				for (my $i=0;$i<=$#header;$i++) {
					if (defined $out[$i]) {
					}else{
						$out[$i] = 0;
					}
				}
				my $out = join("|",@out);
				$out = "|$out|";
				if ($out !~ /type\|num/) {
					$content = "$content$out\n";
					#print "$content\n";
				}
			}
		
		#print $content;<>;
		$report_hash{$type} = $content;


	}

	sub row_merge{
		my $type = shift;
		my $header;
		my @content;
		my $count = 0;
		foreach my $samples (@all_sample) {
			open EACH,"$samples/$report_hash{$type}" || die $!;
			####拿到header
			if ($type eq "mutation_statistic") {
				my $header1 = <EACH>;
				my $header2 = <EACH>;
				$header = "$header1$header2";
			}else{
				$header = <EACH>; ##拿到title
				$header =~ s/^#//g;
			}
			while (<EACH>) {
				chomp;
				my @lines = split("\t");
				$lines[0] = $sample_name[$count];
				my $lines = join("|",@lines);
				push @content,"|$lines|";
			}
			close EACH;
			$count++;
		}
		
		###1. 转化header的格式
			my @header = split(/\t/,$header);
			$header = join("|",@header);
			$header = "|$header|";
			$header =~ s/\n\|/\|\n/g;
			$header =~ s/\|A>C/\|\|A>C/g;
		###2. 制作表格的格式
			my $format = "|";
			if ($type eq "mutation_statistic") {
				@header = split(/\|/,(split(/\n/,$header))[0]);
				pop @header;
			}
			foreach my $run (@header) {
				$format = "$format:-:|";
			}
		###3. 得到内容的格式 
			my $content = join("\n",@content);
		
		$content = "$format\n$header$content\n";
		#print "$content";<>;
		$report_hash{$type} = $content; ###在刚开始的时候,hash是文件的路径，当所有的路径运行完毕之后，统计文件的values就会替换成markdown的表格形式	

}




sub check{
	###获取QC下的ZIP文件###
	my $count = 0;
	my $sample;
	my %hash;
	open LIST,"$list" || die $!;
	while (<LIST>) {
		chomp;
		my @lines = split(/\t/);
		$sample = $lines[0];
		push @all_sample,$lines[0];
		push @sample_name,$lines[1];
		chdir($lines[0]);

			opendir DIR,"1_QC" || die $!;
			my @dir = readdir(DIR);
			close DIR;

			my $R1_fastqc;
			foreach my $dir (@dir) {
				if ($dir =~ /zip$/ && $dir =~ /R1/ && $dir !~ /^QC/) {
					$R1_fastqc = $dir;
				}
			}

			chdir("1_QC/");
			system("unzip -o $R1_fastqc 1>zip.log"); ###解压
			chdir("../");

			$R1_fastqc =~ s/.zip$//g;


			####
			%hash = (  ###给定位置
				"raw_Q20" => "1_QC/raw_statistic_reform.xls",
				"insert_plot" => "4_annot/insertion_len.png",
				"ATCG_per" => "1_QC/$R1_fastqc/Images/per_base_sequence_content.png",
				"base_quality_per" => "1_QC/$R1_fastqc/Images/per_base_quality.png",
				"seq_quality_per" => "1_QC/$R1_fastqc/Images/per_sequence_quality.png",
				"reads_dup" => "1_QC/$R1_fastqc/Images/duplication_levels.png",
				"clean_Q20" => "1_QC/clean_statistic_reform.xls",
				"mapping_statistic" => "4_annot/mapping_statistic.xls",
				"exon_statistic" => "4_annot/exon_statistic_result.xls",
				"snp_indel_statistic" => "4_annot/SNP_INDEL_statistic.xls",
				"snp_base" => "4_annot/SNP_base.png",
				"indel_base" => "4_annot/INDEL_base.png",
				"int_del_plot" => "4_annot/int_del_plot.png",
				"mutation_statistic" => "4_annot/mutation_type.xls",
				"snp_annot_statistic" => "4_annot/snp_statistic.xls",
				"indel_annot_statistic" => "4_annot/indel_statistic.xls",
				"snp_annot_pie" => "4_annot/snp_statistic_pie.png",
				"indel_annot_pie" => "4_annot/indel_statistic_pie.png",
				"snp_pro_statistic" => "4_annot/snp_protein_statistic.xls",
				"indel_pro_statistic" => "4_annot/indel_protein_statistic.xls",
				"sift_statistic" => "4_annot/sift_statistic.xls",
				"hdiv_statistic" => "4_annot/hdiv_statistic.xls",

			);

			foreach my $keys (keys %hash) {
				if (-e "$hash{$keys}") {
				}else{
					print "$lines[0]:$keys do not exist check your file\n";
					$count ++;
				}
			}

		chdir("$pwd_path");
	}
	close LIST;
	if ($count > 0) {
		print "some files are missing\ncheck your program and rerun\n";
		if ($force == 0) {
			exit;
		}
	}else{
		print "check pass!\n";
	}

	my @back = (\%hash,$sample);
	return(\@back);
}


sub makefolder{
	mkdir("5_report",0777);
	mkdir("report_temp",0777);
}

sub delete_windows_enter{
	open IN_SUB,"$list" || die $!;
	open OUT_SUB,">${list}_sed" || die $!;
	while (<IN_SUB>) {
		chomp;
		my $lines = $_;
		$lines =~ s/\r//g;
		print OUT_SUB "$lines\n";
	}
	close IN_SUB;
	close OUT_SUB;

	system("mv ${list}_sed $list");
}